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Abstract: B-mode polarization of the cosmic background radiation is induced from purely 
scalar primordial sources at second order in perturbations of the homogeneous, isotropic 
universe. We calculate the B-mode angular power spectrum Cj^^ sourced by the second- 
order scattering term in the full second-order Boltzmann equations for the polarized radi- 
ation phase-space density, which have recently become available. We find that at / w 200 
the second-order effect is comparable to the first-order effect for a tensor-to-scalar ratio 
of r = 10~®, and to about 2 • 10~^ at / ~ 1000. It is always negligible relative to the 
weak-lensing induced contribution. 
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1. Introduction 

Polarization is expected to play a central role in future studies of the cosmic background 
radiation. The polarization patterns are usually split into a divergence-like E-mode and 
a curl-like B-mode |]^, |2[. B-mode polarization is a powerful diagnostic for primordial 
gravitational waves, i.e. tensor fluctuations of the metric, and thereby constrains inflation 
models directly. While E-mode polarization has already been detected Q and is being 
observed with increasing precision [^, a B-mode signal remains elusive. This, together 
with precise information on the temperature anisotropy spectrum indicates some suppres- 
sion of primordial tensor versus scalar perturbations, since B-modes are not generated by 
scalar perturbations. 
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The absence B-mode polarization when the primordial fluctuations are purely scalar 
holds, however, only in linear perturbation theory. If primordial tensor fluctuations are 
indeed suppressed, B-mode polarization generated from scalar sources in second order may 
constitute an important background to the search for primordial gravitational waves. While 
such an effect would naturally be expected to be relevant at tensor-to-scalar ratios of order 
10~^, which is the size of perturbations in the microwave background, only a full second- 
order calculation can tell whether there are no enhancements. In this paper we compute a 
new second-order effect that contributes to B-mode polarization. 

Several second-order sources of B-mode polarization have already been calculated in 
various approximations. The most important is the weak-lensing effect, reviewed in Ref. Q, 
which converts E-mode polarization to B-mode polarization as the photons travel through 
the inhomogeneous universe |^]. Weak lensing becomes large at small scales, and at large 
values of the perturbation wave-vector k the perturbation series breaks down. The usual 
treatment of weak lensing therefore avoids perturbation theory by considering the small 
deflection angles of the photon trajectories. Another effect that has been estimated is 
B-mode polarization from gravitational time delay Q and from sources proportional to 
second-order vector and tensor metric perturbations, which are themselves generated from 
the product of scalar perturbations [10|. 

Recently the full second-order Boltzmann equations for the cosmological evolution of 



the polarized radiation distribution have become available [11, so that in principle it 



is possible to compute the power spectrum of B-mode polarization in second order gen- 
erated from primordial scalar sources without approximations. While this is numerically 
challenging and beyond the scope of the present paper, we focus here on the novel sources 
of B-mode polarization that appear in the second-order collision term, which have not been 
estimated before. As in previous investigations of the weak-lensing and second-order met- 
ric perturbation effect, we calculate the amount of B-mode polarization from second-order 
scattering sources in isolation, that is, we set the sources of other effects to zero in the 
equations. We then compare the magnitude of the collisional effect to those previously 
known. Investigating the various effects in isolation may be viewed as a first step to solv- 
ing the full equations, and is probably a good approximation, since the different effects 
are related to different periods in the evolution of the universe. In particular, the collision 
sources relevant to the present paper are active mainly during a short period around re- 
combination, and later, during reionization. The reionization contribution will, however, 
not be considered here. 

The outline of the paper is as follows. In Section |2| we define the B-mode angular 
power spectrum when the radiation spectrum is not black-body. We briefly review the 
structure of the second-order polarized Boltzmann equations and discuss the terms that 
define the new collisional sources. A numerical method based on Green functions, which 
is suited to compute the two-point function that relates to the angular power spectrum 
C^^ of B-mode temperature fluctuations, is derived in Section ^. In Section |^ we apply 
this method to calculate the power spectrum and discuss issues related to numerical checks 
and stability. Our results and conclusions are summarized in Section ^. Two appendices 
collect equations and derivations that supplement Sections § and |3[ 
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2. Second-order scattering sources of CBR polarization 



In this section we begin by defining the B-mode angular power spectrum and summarize 
the second-order Boltzmann equations that will be solved subsequently. 

2.1 The B-mode angular power spectrum 

The cosmic background radiation (CBR) is described by the phase-space density matrix 



fab{r], x,q) = 6ab fP (q) + f^^ (v, x,q) + f^J^' {7],x,q) + ... 



(2), 



(2.1) 



which we expand around the unpolarized, homogeneous black-body background distribu- 
tion fj^\q). Expressed in terms of co- moving phase-space momenta q = qn, ff'\q) is 
time-independent and the black-body temperature Tq is the temperature of the CBR to- 
day. The indices a, 6 = it refer to the circular polarization basis and can be exchanged for 
the Stokes parameters X = I, V,E,B hy a linear transformation |11]. We define a matrix 
of fractional temperature perturbations @ab{VjX,q) through 



fabir],x,q) 



exp 



To (1 + 9(7?, a;, q)) 



1 



(2.2) 



where the "one" in square parenthesis and in the argument of the exponential function 
must be interpreted as the 2x2 unit matrix, and we expand Q{ri,x,q) = 0(^)(r/, a;, n) + 
Q^'^\r],x,q) + . . .. The term "temperature perturbation" is slightly misleading, since in 
general fabiv^ q) is not a Bose-Einstein distribution when perturbations are included, and 
hence Q{r],x,q) is not independent of the photon energy q. However, as is well-known, 
the first-order perturbation Q^^\r],x,n) is independent of q, and the spectrum remains 
black-body at first order with a position- and direction-dependent temperature. At second 
order, however, distortions of the black-body spectrum are expected. 

We wish to consider observables independent of q, which may be interpreted as tem- 
perature perturbations. To this end we define the radiation energy density normalized to 
the unperturbed energy density, 



Aab{v,x,n) 



f dqq^fab{v,x,qn) 



Jdqq^fPiq) 



(2.3) 



If the photons were unpolarized and assumed a black-body spectrum with temperature 
T{r],x,n), then A{r],x,n) would be related to the temperature perturbation by 



A{T],x,n) 



T{r],x,n) 



Tn 



il + Q{7],x,n)Y 



(2.4) 



In general we obtain the matrix identities 

A^^\7],x,n) = A&^^\ri,x,n) 



A(2) (r/, x,n) = ie^ ' (77, x,n) + 6 Q^^^ (r?, x, n) 



(2.5) 
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among the perturbation coefficients, where 

fdqq^ff\q) \Q(^){r^,x,qn) + lqM^){rj,x,qn) 
e^ ' {t], X, n) = 1 ' 



jdqq^ff\q) 



(2.6) 



The quantity 0^^'' obviously equals 0^^^ when the latter is q independent. Otherwise it 
represents the temperature perturbation of a black-body distribution with the same energy 
density as the actual radiation distribution f{r],x,q). 

We define the multipole expansion coefficients of the fractional temperature perturba- 
tion here and today as usual by 



dn{n)YiZ{n 



e(i)(7?o,a3o,n) + e^'^(7?o,a3o,n) + ... , (2.7) 



noting the presence of the equivalent black-body temperature perturbation and the 2x2 
matrix nature of aim- The spin of the spin- weighted spherical harmonic, Y^^(n), is s = 

for the diagonal ++, components, s = +2 for ab = + — , and s = —2 for — h- In the 

circular polarization basis B-mode polarization is given by i/2 times the difference of the 
-\ — and — h components of the 2x2 density matrix, hence we define 



0-B,lm — 2 ('^H — ~ 0- — l-,lm)- 



(2.8) 



The absence of B-mode polarization in first order implies that the H — and — -|- components 



fab^ ®ab ^^^^ [®*'^^]ab equal, SO 



^^mm = 0- At second order 



'■BAm 



dn (n) (n) {m,xo,n)- y-J* (n) {m ,^o,n) 



d^k 
(2^ 



Jkxo 



^f^^BMiVo,k), 



(2. 



where in the last line we introduced the Fourier modes and multipole coefficients of 
A2^(r?o,a;o,n). 

Our aim is to compute the B-mode angular power spectrum C^^ given by the statis- 
tical average 

{0'B,lm^*B,l'm') = ^U'^mm' (2-10) 

when the perturbations of the FRW background at first order are purely scalar. Using 
Eq. ( |2.9| ) this average can be expressed in terms of 

(2.11) 



in the form 

{o-B,lm^*B,l'm' 



J_ .y 4' 47r 

16 ^ ' ypTTWTl) J (27r 



(2-12) 



The power spectrum of As defined here depends on the direction k of the mode vector, 
since it refers to the fixed coordinate system of an observer and helicity m defined with 
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respect to the three-axis of this system. When the heUcity axis is chosen to be k we obtain 
the simpler expression 

(Ag,L(%>fc)Agiw(r/o,fc'))|fcaxis = W (2^)^5(3) (fc - k')P,^l^{k) (2.13) 

where the form of the right hand side of the equation follows from statistical isotropy and 
homogeneity and explicitly from the results of Section ^. In particular, the power spectrum 
-^z^m(^) depends only on the magnitude k = \k\. It is this power spectrum that will be 
computed later on for I' = I. The two power spectra are related through the transformation 
( |A.5| ) of the B-mode multipoles under rotations, which gives 

Pw%nAk,k) = ^(^^ ^ Y.'^rrtik)Yf^\k)P^,-^{ky (2.14) 

Plugging this into Eq. ( |2.12| ) and using the orthogonality of the spin-weighted spherical 
harmonics, we obtain the familiar result 

(2) 

The sum is restricted to values |m| < 2, since the equations for /\'-^ i^{k) in the frame where 
k coincides with the helicity axis do not contain source terms when |m| > 2. The sum does 
not include m = since there are no scalar mode contributions to B-mode polarization. 

2.2 Second-order equations 

In Ref. |ll| we derived the second-order Boltzmann equations for the photon phase-space 
densities fx^imiVi^iO)- We now explain the structure of these equations and discuss the 
approximations we apply in this paper to isolate the second-order scattering sources. We 
are interested in the energy integrated distribution functions defined by 



,(2) 



in particular in A^ ^^(7/, fc). The equation for second-order B-mode polarization when 
there are no first-order vector and tensor perturbations, hence A^l (??, k) = 0, is given by 



mim 



+ (AW-Z?«)(fcOAi5,L^(fe2)i^° 
+ {^k^r^ (^(^) - D('^) (fei) + 4.fcr-U(i)(fci)) Ag^^(fc2X^ 

- X^S-]^^^^ (^S™^ - -^'<2miyk2)D^mlm}- (2.17) 
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Some comments on notation are in order. The equations above and below are given in 
conformal Newtonian gauge for phase-space densities defined in an inertial frame locahy 
at rest and ahgned with respect to the general coordinate system. A^^^ and D^^^ denote 
the first-order scalar metric perturbations, We,[m] the velocity field of the baryon-electron 
fluid, and \k\ is proportional to the collision rate of Thomson scattering. The coupling 
coefficients Dm[m etc. are summarized in Eq. ( |A.1C| ) and Ref. [11|. Products of first- 



order perturbations involve two different mode vectors ki and k2 = k — ki and it is 
understood that a convolution integral over ki is performed. Similarly, for given m we 
imply 7712 = 771 — 77ii and a sum over 7711. Indices in square brackets refer to helicity 
rather than Cartesian vector components with k^'^^ = ik^. We refer to Ref. |11] for further 
notational details. 

Eq. ( |2.17D and the corresponding equations for the radiation intensity and E-mode 
polarization contain a number of different effects: 

• In addition to the conformal time derivative the first line contains the effect of free 
streaming of radiation perturbations in the unperturbed background. This term 
converts vector and tensor E-mode perturbations into the corresponding B-mode 
perturbations and further excites higher multipoles from the I = 0,1,2 ones after 
radiation ceases to be tightly coupled to the baryons. These terms must obviously 
be included in a calculation of the angular power spectrum observed today. 

• The remaining two lines before the equality sign are the weak lensing and time delay 
contributions, which represent the effect of space-time inhomogeneity on the photon 
perturbations as they travel through the universe. This converts E-mode polarization 
into B-mode polarization even in the absence of vector and tensor perturbations and 
is known to generate significant B-mode polarization at large I We drop these 
known contributions, since we are interested in the effect of the collision source terms. 

(2) 

• The equation for Aj {^{k) contains second-order metric source terms involving prod- 
ucts of first-order perturbations such as A^'^\ki)D^^\k2) and second-order perturba- 
tions A'^'^\k), B^T^{k) and E^T^{k). The amount of B-mode polarization generated 
from some of these terms has been estimated in Ref. Iin]. Specifically, the vector and 
tensor perturbations ByJ^-^{k) and EyJ^-^{k) induced at second-order from purely scalar 
perturbations have been evaluated and inserted into the first-order equations for free 
streaming to obtain the amount of B-mode polarization today. Again, since we are 
interested in the effect of the new collision source terms, we drop these terms. 



The right-hand side of Eq. ( 2.17| ) is the second-order collision term. The first term 



in curly brackets is a universal relaxation term for all multipoles that in the absence 
of source terms drives the phase-space distribution to its equilibrium form and hence 
damps polarization. This term is already present in first order. In addition we find 
new scattering sources of B-mode polarization from the coupling of first-order scalar 
intensity and E-mode polarization to perturbations in the electron-baryon bulk ve- 
locity. These terms and corresponding new terms in the equations for Aj l^{k) 
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and A^ ^,^^^(fc) arc kept in our numerical calculation. They constitute the scattering 
sources, whose effect is computed here for the first time. We thus include the com- 
plete second-order collision term with one exception. We neglect terms of the form 
\k\ [Sxe/xe]^^^ x first order perturbations, which arise from the perturbed ionization 
history. 

We thus solve the following system of second-order equations: 

± 

+ ^[U] (fci) (a?L. - V6Ag),„ J (fc2)C-f„ } (2.18) 



± 

- '^'3^<[U] (fci) (Ai;L. - v^Ag),^ J (fc2)I^-f„ } (2.19) 
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,(2) (U\ 1 .,(1) It. \a(1) (U \ T^O,l 



- 6,2^v^Uk,) (Ai^L, - ^^A« „J(fc,)Z):^?„ } (2.20) 

+ <[L](^i)ASo(fc2)}. (2.21) 
with i? = 3/>b/(4/0'y). The equation for 'ygr^i(fc) must be included, since this quantity 



appears in the colhsion term of the intensity perturbation equation ( 2.1^ ), and we apphed 



the same approximations to this equation as for the other three. That is, we neglect 

(2) 

the second-order metric perturbation -Bj^j , and for consistency also products of first-order 

(2) 

metric perturbations (since is itself sourced by a product of such terms). The complete 

('2') ('2') (1^ 

equation for and the baryon density perturbation 61 = [dpb /pb] and the precise 

definition of these quantities is provided in Appendix^. Eqs. ( p. 181) to (|2.21| ) form a closed 
system together with the equations for the first-order perturbations also summarized in the 
Appendix, which we solve in the following. We note that the approximations made are 
systematic in the sense that in the absence of collision terms {\k\ — )• 0) all second-order 
perturbations vanish; hence we exclusively focus on collisional effects as intended. 

3. Solving the second-order equations 

For convenience, we introduce a compact notation summarizing the photon equations 
Eq. (|2T8D to Eq. ( ^ ): 



Ai^) + kCn^^'l^ = -|/i|(Ai2) - ^„^A(^) - Sn). (3.1) 

The indices are multi-indices, n = {Xn,lnmn), where X = I, E, B, Ve distinguishes between 
the photon modes and the baryon velocity; I, m are the multipole indices. Repeated multi- 
indices are summed over X, I, and m; k has been aligned with the three-direction, k = ke^. 
The dependence on k of A„, Sn, and the Green functions Gnm introduced below will not 
be made explicit. Note that Eq. ( |3.lD does not represent the electron equation Eq. ( 2.21J ), 
i.e., it is only valid for n ^ (ve, Im). 

The matrix Cnm describes the free-streaming coupling of each photon multipole mo- 
ment to its neighbours with / it 1, which leads to a gradual excitation of the initially small 
large-/ moments. The conformal time it takes for an excitation to propagate from multipole 
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moment / to moment I ziz Al is 77 ~ Al/k. Free streaming also accounts for the conver- 
sion between E- and B-mode polarization for vector and tensor modes, which is the only 

possible source for B-mode polarization at first order. 

(2) 

The first part of the scattering term, — |K|An , is responsible for the tight-coupling 

suppression. At early times the scattering rate \k\ is large, so that a non-vanishing moment 
(2) 

An induces a large gradient driving it to zero as long as the first part of the scattering 

(2) 

term is not cancelled by the remaining term <^nm^rn + Sn- Since the coupling ^nm vanishes 
for high multipole moments, only the monopole and dipole are not suppressed by the tight 
coupling of photons and baryons. At second order the quadrupole A/ 2m is also present 



due to a cancellation with Sn but there is no polarization in tight coupling 

The source 5„ contains the convolutions of first-order perturbations. Other than in 
the first-order equations, where there are only source terms for low multipole moments, in 
the second order source Sn there are contributions for all multipoles. Moreover, there is a 
source term for B-mode polarization which is induced by first-order quadrupoles, including 



the large intensity quadrupole, cf. Eq. ( p. 20 ) 



Sbm = ••• -fe^i;l'L,](^i)^i!L,(fc2)I?^?m--- • (3.2) 

Such a direct source term for B polarization does not exist at first order, where B-mode 
polarization is only generated from free streaming of vector and tensor E-mode polarization. 

3.1 Line-of-sight integration 

At first order the line-of-sight integration is used to great success ^ . It is the solution 
of 

A„ + kCnrn^rn — ~\^^\^n + Pni (3-3) 

where Cnm are the free-streaming coefficients introduced in Eq. (^]^) and pn includes all 
source terms. The solution to this equation is given by 

An(r?o) = r dve-'^'-'^^jnmiMm - v))pm{v), (3.4) 



JQ 

where the functions jnm are combinations of spherical Bessel functions and Clebsch-Gordan 
coefficients, 

jn^x) = Y,^'--'^-'- t/'^f 1 ^ 3h{x)Hx„xJln -h- Im) 




(3.5) 



The matrix Hxx' is responsible for the mixing between E and B polarization in free 
streaming and is given by Hxx'i^^^^) = ^xx' and Hxx'{odd) = 5xiSx'i + i^XE^X'B — 
iSxsSx'E Fx = when X = I and Fx = —2 when X = E, B. Clebsch-Gordan 



coefficients are denoted by big round brackets. A derivation of this solution is given in 
Appendix |^. The functions jnm are always real since the product can 
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be imaginary only for odd In — h — Im and X„ = Xm = I, but in that case the second 
Clebsch-Gordan coefficient vanishes. 

To employ this solution at second order we only need to use the appropriate source 
terms. In Eq. (|3.1|) we identify pn = |k|(S'„ + <;nm^rn) and we obtain the integral solution 



rvo 

^^nHm) = / dv |K(7?)|e-'^('')j„„(A:(r?o - + )(^)- (3-6) 



Despite the second-order term in the integrand, this is a big step towards the computation 
of second-order quantities today: the coefficient vanishes for multipole moments Ip > 2 

(2) 

so that we are left with the task of finding solutions for Ap (rj) with I < 2. In addition, 
the visibility function |K;|e~'' is non-zero only around recombination and thus it is sufficient 
to compute these solutions for early times rj < 500Mpc/c. (We quote conformal time in 
units of Mpc/c, but set the speed of light c = 1 in equations.) 

3.2 Solution using Green functions 

The Boltzmann equations have to be solved with stochastic initial conditions resulting 
from inflation. At first order, the linearity of the equations allows to write the solution as 
product of transfer functions and the primordial fluctuations. The problem reduces to the 
task of computing the non-stochastic transfer functions. This straightforward separation 
is no longer possible at second order due to the quadratic convolution terms. However, we 
can achieve a similar separation by using Green functions. 

The second-order photon equation Eq. (|3.1|) is linear in the second-order quantities, 



A(2)(,;)=^_(^)A(^)(r,) + a„(r?), (3.7) 

where Anm = |^| (?nm — Snm) — kCnm and an = \K\Sn- The same applies to the electron- 
velocity equation with Anm and o"„ according to Eq. ( p. 21 ). The solution of such a linear 
differential equation can be written in terms of Green functions, 

rv 

^n\v) =Gnrn{r], Vini) ^'fn'i'Hmi) + drj' Gnm.{V,V')<^m{v)-: (3-8) 



where the Green function GnmiViV') = Gnm{v->v')^{v ~ v') satisfies 

dr^GnmiV^v') = ^np(??)^pm (?/, ?/') + Snm^iv " v')- (3-9) 

In this differential equation there are no stochastic quantities, thus we can compute the 
Green function GnmiViV') by solving the equation 

dr^GnmiV, v') = Anp{r])Gpm{v, V') (3-10) 

for rj > rj' with initial condition Gnmiv' iV') — ^nm numerically using standard methods. 
Note that with k = ke^ aligned with the three direction and ik^^^ = —k, the matrix Anp{r]) 
is real and therefore the Green functions are real as well. 

In principle, Green functions could be used to calculate all two-point functions for 
any r] until today, but generating the Green functions for late times and large multipoles 
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numerically is very time-consuming. A much better performance can be achieved by com- 
bining this ansatz with the line-of-sight integration. Using the Green function method we 
compute Ap'^^(77) for / < 2 and for early times, and then use Eq. ( |3.6[ ) (with lower limit 
r] = replaced by r/ini, since the visibility function is negligibly small at very early times) 
to obtain the late-time evolution of the higher multipole moments: 

rvo r 
^^nHvo) = / dr]\k{r])\e-^^'^^jnm{HvO-V)) <impGpq{ll,nmi)^f\r]ini) 

J Vini 

r' r I /I 

+ / dr]' (6rnq5{il - ri') + \k{r]')\qrnpGpq{ri,ri')) Sq{'q') . (3.11) 

J Vini 

Note that besides Ag (r/ini) only the source Sq{r]') depends on the stochastic initial con- 
ditions whereas the evolution is given by the non-stochastic functions jnm and Gpg. This 
separation is convenient for computing second-order correlation functions as described in 
the following subsections. 



3.3 Second-order pow^er spectra 



The result Eq. ( p.ll ) simplifies further if the initial second-order quantities Ag^^(?7ini) van- 



ish. As we will discuss in section this is the case for the computation of the non- 
scalar modes with m = ibl,ib2 and only these are needed in Eq. (2.15) to compute the 



B-mode angular power-spectrum. Consequently, we evaluate the two-point function as- 
suming Aq (ryini) = and obtain 

(A(f)(fc,ryo)Ag)*(fc',r?o)) 

d7? / dr?' \kirj)\\kirj')\ e'^^^^^'^^^ ^jnmiHvo - r?))jn'm'(M% - v)) 



Vini -J Vi 



rv rv 

X I dr]i dr][{6mqS{r] - Tji) + \k{r]i)\<^rnpGpq{r],r]i)) 



Vini J Vi 



X [dm'g'Siv' -V'l) + \kiv'l)U'p'Gp,q,{v',v'l)){Sqik,m)S*q>ik',v'l)), (3-12) 

where we anticipated Eq. ( p.l6| ) below that sets fc = fc' in the integrand. We also imply 
that the helicity axis is the direction of fc or, equivalently, that k is aligned with the three- 
axis of the coordinate system. Comparison with Eq. (|2.13|) provides the power spectrum 
of A^ required to compute the B-mode angular power spectrum (|2.15| ). The source term 



correlation function {Sq{k,r]i)S*,{k' ,r][)) can be calculated from first-order quantities only. 
These are given by the primordial potential $(fc) = A^^\k,r]i^[), multiplied by a transfer 
function T. Using k2 = k — ki we can write 

^i^rA«(fci,r/)AW(fc2,r/) 

^ KPJ T^i) (fci, 7?)r(i) {k2, vMkiMk2), (3.13) 
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where the constants Kn"^ are determined by the expressions for the second-order colhsion 
term, Eqs. (|2.18| ) to ( 2.21| ). Thus, to compute the expectation value on the right-hand side 
of Eq. ( |3.12| ) we need to evaluate a four-point function of the primordial potential 

= mk,)^k2)){^*{k[)^*{k',)) + mki)<^*{k[))mk2W{k',)) 

+ {<^{k^)^*{k',))mk2W{k[)) 

= {2Trf(^P^{ki)P^{k[)6'^^\k)6^^^k') 

+P^{h)P^{k2)\s^^\ki - k[)5^^\k2 - fc^) + - fc^)5(3)(fc2 - k[)]), (3.14) 



where P^ is the primordial power spectrum, (^>(fc)$*(fc')) = {2TTf6^^\k-k')P<s>{k), and 
Gaussian statistics of the primordial potential $(fc) was assumed. For fc 7^ or fc' 7^ this 
results in 



mki)^{k2W{k[)^*{k'2)) = {2T,fP^{ki)P^{k2)5^^'^ (k - k') [5(3) {ki -k[) + 5(3) (fci - k'2 

(3.15) 



and the complete expectation value on the right-hand side of Eq. (3.12) reads 



(54fc,r?)5:,(fc',7?0) = (27r)3<5(3)(fc-fc') 



(fki_ 
(2^ 



P^{ki)P^{k2)KPiKP,'^ 



T«(fci,r?)r«(fc2,r?)r«*(fci,r/')Ti,^^*(fc2,V) 



- r(i) (fci, 7?)r(i) {k2, r?)T«* {k2,v')T^}^* {ki , v' 



(3.16) 



For a given primordial spectrum, we can now compute the second-order power spec- 
trum using Eq. ( 3.12| ) and Eq. (|3.16|) . 



3.4 Non-Gaussianity 

Analogously, the method described above can be applied to study non-Gaussianity. At first 
order, the bispectrum (three-point function) is always proportional to the primordial bis- 
pectrum and therefore zero if the primordial perturbations are Gaussian. This is different, 
if non-linear effects are taken into account. The leading contributions to the bispectrum 
are terms combining one second-order perturbation with two first-order perturbations. 



(A„A^Ap) = (Ai2)AWAW) + sym. 



(3.17) 



We replace the first-order quantities by their transfer functions and the primordial potential 
and calculate the second-order quantity using the combination of line-of-sight integration 
and Green functions, Eq. ( p. 11 ). Allowing for non-vanishing initial values An\k,r]ini) we 



- 12 - 



find, for k aligned with the three-axis of the coordinate system, 
(A(2)(fc,r/o)AW(fc',r?o)AW(fc",7?o)) = r«(fc', 



X / dr]\k\e''^jnmikir]o - r])) 

'Vini 

J Vini 



. (3.18) 



Using Eq. ( ^.13 ) and contracting the resulting four-point function similar to Eq. ( 3.14 ) we 
can write the second expectation value as 



{Sn{k,r,')ct>{k')cl){k")) = {2Trf6^''Hk + k' + k")P^{k')P^{k")KPi 

X rji) i-k', v')tW i-k", v') + tW v)r(i) i-k', rj'] 



(3.19) 



(2) 

For an initial value Aq (fc,?/ini) which is quadratic in the primordial potential ^, also the 
first expectation value can be contracted and written in terms of the primordial power- 
spectrum. A detailed study of non-Gaussianity from second-order effects employing these 
equations will be presented in a follow-up article. 



4. Numerical evaluation 



In this section we outline the steps required to compute the two-point function ( |3.12| ) 
and the C^^ angular power spectrum ( 2.15| ) numerically. First we need to compute the 



first-order transfer functions, which is discussed in Section 4.1, and the Green functions. 



as discussed in Section 4.2. In Section 4.5 we obtain the source terms and combine both 



results to perform the line-of-sight integral together with the wave- vector convolutions. 

For all intermediate results presented in this section we use a ACDM model with 
massless neutrinos and the following parameters: Tcbr = 2.726 K, Hq = 70km/(sMpc), 
and r^cDM = 0.245, ilbaryon = 0.045, $7a = 0.71. We assume a scale-invariant primordial 
power-spectrum P$(A;) = 27r^(9/25)A^/A;'^ with amplitude A^ = 2.41 x 10~^. To compute 
the recombination history we employ Recfast [^] with a helium mass-fraction Yp = 0.24. 
Our calculation does not include late effects such as the late integrated Sachs- Wolfe effect 
(ISW) and reionization. Scattering effects and the early ISW will be considered up to r/ = 
500Mpc/c and neglected for later times where we use the free-streaming approximation. 
Recombination occurs around t^i-cc = 286.7 Mpc/c, defined as the median of the visibility 
function. 

4.1 First-order solutions 

The first-order transfer functions can be computed by solving the first-order Boltzmann 
equations with standard methods for solving ordinary differential equations (e.g. as imple- 
mented in the GNU Scientific Library [|l7|). Alternatively one can resort to programs like 

or Cmbeasy P^]. The following results are based on our own 



Cmbfast 0, Camb m 
code that has been compared to Camb. 
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A crucial point is that we need only the low multipoles at early times where the 
numerical calculation is straightforward. The restriction to early times is due to the factor 
|k| in all collision term sources which is negligible after recombination. The reason why only 
low multipoles are needed at early times is that higher moments are only slowly generated 
by free streaming. Thus, one can cut the Boltzmann hierarchy at some /cut and neglect 
all higher multipole moments. Alternatively, to save CPU time, one can apply the cut 
at a much lower multipole where the corresponding moment cannot be neglected, but is 



replaced by a closing relation. This option will be discussed in detail in Section |4.4 . 

Neutrinos start free streaming much earlier than photons because they are not tightly 
coupled to baryons by Thomson scattering and thus start generating higher multipoles 
at early times. This has to be taken into account by choosing a cut at larger I for the 
neutrino hierarchy than for the photon hierarchy. Typically we take more than 50 neutrino 
multipoles into account which yields accurate results until matter domination. 

4.2 Initial conditions 

We assume that the primordial first-order perturbations are adiabatic and begin the evo- 
lution of the transfer functions at ajni = 10^^ (corresponding to rji^i = 0.464 Mpc/c) deep 
in the radiation era with standard adiabatic initial conditions for the scalar perturbations: 



5p^ + p, 



^ZUV^) = TSUm = \T'»Jk) = '-Tjiljk) = (4,1) 

For the neutrinos we further include an initial quadrupole T^^^^^^{k) = 2k'^/{3H^), which 
however is very small at ?7ini, when for all k of interest k/Hc <C 1. 

Setting the initial conditions for the second-order perturbation variables is simplified 
by the fact that by Eq. ( p.l5| ) we only need to compute the vector and tensor perturbations 
m = ±1, ±2 when working in the frame where the mode- vector k is aligned with the helicity 
axis (usually the three-direction). While we may adopt the convention that at some initial 
time A^"^^ vanishes, since any non-zero value can be absorbed into a small change of A^^\ 
this cannot generally be done for all perturbation variables. For instance, deep in the 
radiation era, when k/Hc <^ 1, the total energy density perturbation is given by 

^ = _2^(2)+4i^^[^(i)]2. (4.2) 
P 

However, the equations we need to solve for m = ibl,ib2 do not contain the second-order 
energy density perturbations, that is the monopoles of the multipole decomposition, and 
hence we do not have to determine their initial conditions. 

(2) 

We do need the m = ±1 components of the second-order electron velocity v , , and the 
radiation intensity dipole. Their initial values are related to the metric vector perturbation 

(2) 

B}{ . Since, as discussed above, we do not consider the second-order metric perturbations 
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in this paper, the initial value of the second-order velocities and dipoles is consistently 
set to zero. In the tightly coupled radiation era there exists a non-vanishing radiation 

(2) 

quadrupole Aj which acquires a non-zero initial condition given in terms of the square 
of the first-order electron velocity. Being of second order and suppressed by (k/Hc)'^, this 
can be safely neglected. We show in Section 4^ that collisions quickly drive the quadrupole 
to its tight-coupling value when it is zero initially. Finally, we note that the second-order 
neutrino perturbations do not appear in our equations, so we do not have to set their initial 
values. 

To sum up, we may solve the second-order equations with all initial values of second- 
order variables set to zero. 



4.3 Green functions 



The differential equations ( 3.10 ) for the Green functions are far less complicated than the 
original second-order equations. In particular, the equations are no longer stochastic. Their 
structure is identical to the one of the first-order Boltzmann equations. 

The relevant Green functions Gnmiv^ v') can be classified by several criteria. All Green 
functions except those for monopole, dipole, quadrupole and electron velocity are strongly 
tight-coupling suppressed.^ This leads to a suppression of these Green functions if rj' is 
located before recombination (except for r] close to t/'), and allows us to restrict the ij' 
integration to a period directly around recombination, as after recombination the source 
terms vanish with |k| and before recombination the Green function is small, even after 
multiplication with the large scattering rate. The effect resembles that of the visibility 
function l^le"*^ in the line-of-sight integration which is peaked around recombination. 

Green functions which are not suppressed during tight-coupling may lead to numerical 
difficulties because they in turn do not suppress the source an = |k|S'„ in Eq. ( |3.8|) . At 
early times, the scattering rate \k\ is huge and, multiplied with the numerical error of 
Sn, yields a large absolute error of (T„ and the corresponding second-order quantities. As 
a consequence, linear combinations of these second-order quantities can obtain a large 
relative error, if there exist cancellations. For the electron velocity and radiation dipole 
Green functions such problems can be avoided by exploiting the close relation of dipole 
and electron-velocity sources. The electron-velocity sources can be split into two parts, 
one which is - apart from a prefactor ~ identical to the photon-intensity dipole source and 
a remaining part Sy^^im which cancels the factor of \k\ by the tight-coupling suppression 
1/1 k| of the combination of first-order perturbations multiplying it. Writing 

Sve,lm — Sj^im + S-u^^imi (4-3) 



where 

^ i_ = 

4R ^'"^^ V ""e.^m 



S..,in. = i Sl'\k,) Uv['Uk,) - AfUk,) ) , (4.4) 



""^This does not apply to the corresponding calculation of second-order neutrino perturbations. 
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Figure 1: In dashed (red) the Green function G'(£;^2i),(/,ii)(500 Mpc/c, ry') for fe = ke^, 
k — 0.02 Mpc^^. Since the dipole is not tight-couphng suppressed, the Green function does 
not vanish before recombination. After recombination the Green function vanishes because As, 21 
and are only coupled in the scattering term. The solid (black) line is the combined dipole 

and electron- velocity Green function G(£;^2i),(/,ii) (500 Mpc/c, 77') — 47jG(£;^2i),(t)e,ii)(500Mpc/c, ??') 
which vanishes at early times. 



the first part on the right-hand side of Eq. ( |4.3| ) can be combined with the photon source 
in Eq. ( |3.8| ) to obtain 



^nH??) = / (G'n,(/,lm)(r?,V)'S'/,lm(??') + G'„,(i,e,lm)('?,??')'S'i)e,lm(??')) 

•^»?ini 

+ r dv'\f^{v')\Gnp{^,^')Sp{v') 

all other p ''i^i 

+ / dr]' \k{ri')\Gn,(v,,im){v,v')Sv,,im{v') 

■J Vini 

+ Y r dv'\k(v')\Gnp{ri,v')Sp{v'). (4.5) 

all other p ''i^i 

Thus, only the combination Gn^{i^im) ~ j]^Gn^(ve,im) is multiphed with the source Sj^im 
which is large at early times, but this combination of Green functions does vanish in the 
tight-couphng regime. In Fig. |^ the combination is plotted illustrating the suppression 
during tight-coupling (solid line). 

The remaining Green functions which do not vanish in tight-coupling after multiplica- 
tion with \ k\ include: the Green function acting on the photon monopole, the one acting on 
the remaining part of the electron- velocity source S^^^im and the Green functions acting on 
the quadrupoles. The first does not enter our computation since the photon monopole does 
not couple to polarization, the second is suppressed by Sy^^im which itself is tight-coupling 
suppressed as explained above. Finally, the Green functions acting on the quadrupoles are 




Figure 2: Green functions G(£;^22),(£;,22)(500Mpc/c, ry'), 10 x G(/^22),(£;,22) (500 Mpc/c, 77') and 
G'(/,2i),(£;,2i)(500 Mpc/c, 77') (solid/black, dashed/black and dotted/orange, respectively) for k = 
fees, k — 0.03 Mpc""'^, multiplied with |k(?7')|. The multipole moments AE,2m and Ai^2m are only 
coupled in the scattering term. As expected, the Green function connecting these moments van- 
ishes faster than the Green function connecting moments coupled in the free-streaming term. The 
Green function acting on quadrupoles for \m\ = 1 is only suppressed by \k\ and thus the plotted 
combination is constant at early times. The visibility function divided by two is shown in grey for 
comparison. 



only suppressed by one power of |k| for \m\ = 1, and do not vanish when multiplied by 
Hence, while for \m\ = 2 the t]' integration in Eq. ( |3.8| ) can be restricted to times around 
recombination as there are no contributions at early times, the time integration cannot 
be cut for \m\ = 1. This difference between vector and tensor perturbations is related to 
the persistence of the dipoles during tight-coupling. While during tight-coupling no po- 
larization is generated, the short-lived excitation of a quadrupole will modify the dipoles. 
This dipole is then converted into a quadrupole during recombination which couples to 
polarization. This indirect coupling of early vector sources can induce polarization while 
early tensor sources have no influence, as there is no unsuppressed moment. 

After recombination, the characteristics of a Green function depend on the coupling 
between the Green function's two multipoles. If both moments are only coupled in the 
scattering term, the Green function quickly vanishes after recombination. If they are 
coupled via free streaming, it oscillates, since free streaming converts neighboring moments 
into each other. In combination with the factor Eq. ( p. 11 ) both types of Green 



functions decay after recombination. However, the former decays with an additional factor 
of I 1 5 clS depicted in Fig. y. 

4.4 Closing relations 

Like the first-order Boltzmann hierarchy, the equations for the Green functions need to be 
truncated at some multipole moment /cut- The straightforward approach is to set higher 
moments with I > Zcut to zero. This is a good approximation at early times because all 



^The Green functions for |m| = 2 are exponentially suppressed as there is no coupling to the unsuppressed 
dipole. 
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higher moments are initially zero. But once the moment Zcut is excited by free streaming, the 
solutions for the Green functions and first-order quantities, respectively, become inaccurate, 
since in the full hierarchy a non-zero moment at Icut excites the moment /cut + 1 which then 
feeds back into the lower moments. However, the error will first become noticeable only in 
the moments close to /cut- It then takes an equal amount of time to carry the error back 
to the lowest multipoles as it has taken to excite /cut in the first place. Hence, one can still 
trust the results for the lowest multipoles some time after the highest multipole has been 
excited. 

For accurate results one has to cut at some sufficiently large /cut such that the lowest 
moments (which we are interested in) are not influenced by the truncation. To compute 
the photon multipole / at time r] and wave number k, the rule of thumb is to cut at 
/cut ~ I + k (rj — rjj-cc)/2. Thus, if A; or is large, many multipoles have to be taken 
into account. This can be avoided by using appropriate closing relations at a comparably 
small / |14|. 



Approximating the visibility function by a delta function at r/rec allows us to derive an 
analytic relation between the highest multipoles we consider. In this approximation the 
scattering rate is infinitely large before recombination so that all multipoles except for the 
monopole and dipole are zero due to tight-coupling suppression. After recombination, the 
scattering rate and therefore the sources are zero and we can employ the free-streaming 
solution Eq. ( B.lOj ) to obtain 



«2=0,lX',/i ^ 



/ /l /2 W / h h 

m m j \ Fx' Fx' 



H*XX'i^ -h- l2)Ax',l,miVrec)- (4.6) 



Starting from this equation we perform the following steps to derive the closing relations, 
exemplified here for the intensity multipoles A/ iq and /cut = 9-'^ 



We write down Eq. ( |4.6| ) explicitly for the two highest multipole moments below /cut; 

^X,(/cut-l)m and AxXkut-2)m- 

A/,8o(??) = (8j7 - Sjg) A/,io(r/rec) + 17j8A/,oo(^rec) 

A7,7o(??) = (7j7 - 8h) A/,io(??rec) + 15j8A/,oo(?/rcc) (4.7) 

For brevity we leave out the argument of the Bessel functions which is always k{r] — 
r^rec)- For m / 0, the monopole (and even the dipole for m = ±2) does not exist, 
and one uses instead the two lowest non-vanishing multipoles. For polarization, one 
uses the E- and B-mode quadrupole. 

We solve both equations to write the monopole and dipole at recombination in terms 
of the multipoles Ax,(/^m-i)m and ^x,{icut-2)m- In case of polarization {X = E,B), 



^The m — calculation is not required to compute B-mode polarization, but the expressions are less 
complicated and the method can be extended straightforwardly to other values of m. 
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we have to relate AE^(i^^^-i)m and to AE,2m and AB,2m instead, due to 

the mixing of E modes and B modes during free streaming. 



In Eq. ( [4.6D for Ax^i^^^m we can now replace the monopole and dipole on the right- 
hand side with the result from the previous step to obtain the closing relation. In 
our example we find 

. .^ 19i9 (-7j6 + 8js) + 15j7 (9j8 - lOjio) ^ 
17j8 (-7j6 + Sjs) + 15j7 (8j7 - 9j9) 

, 19j9 (8j7 - 9j9) + 17j8 (-9j8 + lOjio) . . ^ i,^. 



17j8 (-7^6 + 8J8) + 15j7 (8i7 - 9J9) 



• Finally, we use this result to replace A^.i^mm in the Boltzmann hierarchy and obtain 
a closed system of differential equations for all moments up to Ax^[i^^^-i)m which is 
independent of higher moments. 

The quality of this approximation depends on the width of the visibility function. While 
the photon evolution is dominated by scattering effects, neutrino modes are exclusively 
sourced by metric terms. Consequently, the closing relations described here cannot be 
applied to neutrinos, since they rely on the sharply peaked visibility function. However, 
we only have to consider neutrino perturbations at first order, so that all Green functions 
can be computed using the closing relations. 

In our calculation of the Green functions we cut the hierarchy at lent = 9, obtaining 
the closing relation 

A/,9m(??) = hl^8m{k{r] - r/rec))A/,8m(??) + hljrn{k{r] - r/rec)) A/jm (r/) . (4.9) 

In general the functions hx,im are combinations of spherical Bessel functions as in the 
171 = example, Eq. ( [4.81) . For m 7^ the functions become more complicated, but as long 
as the argument A;(7/ — ??rcc) is small, they can be approximated by polynomials. This is not 
possible for larger arguments where the functions oscillate, see Fig. ^. Only for very large 
arguments the oscillations can be neglected due to damping and a simple approximation 
can be used again. 

Figure Q shows that the closing relations can be used to significantly reduce the error 
compared to a simple truncation — at almost no additional CPU time. 

4.5 Source terms and integration 

To compute the angular power spectrum Ci at second order, we have to perform eight 
integrations: two time integrals from the line-of-sight solution, another two time integrals 
from the Green function ansatz, three integrations from the convolution over ki and finally 
the k integral in Eq. ( 2.15| ). The integral over the angular coordinate (pi of the wave vector 



ki can be performed analytically. The remaining seven integrals are computed in one single 
Monte Carlo integration. The wave- vector integrations need to be cut at some femax- The 
larger the choice of fcmax; the more multipoles have to be considered when computing first- 
order results and Green functions, increasing the demand for CPU time. Fortunately it 
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Figure 3: The functions /1/.72 (solid/black;) and /i/,82 (daslied/red) wliicli appear in tlie closing 



relation (4.9). 
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Figure 4 

the closing relations (sohd) for A^gj, (black) and Ajjq (orange); k — 0.1 Mpc , icut 
improvement achieved by using closing relations is even better for smaller wave numbers 



turns out that the integrand decays quickly for ki or k2 larger than k. The ki dependence 
of the integrand is illustrated in Fig. ^ depicting the suppression for large values of ki. 
This effect is due to the two primordial power spectra in the expectation value of the 
source terms, Eq. ( |3.16 ). Since k2 = \k — ki\, ioi large ki the power spectra suppress the 
integrand with 0{k^^). We find that choosing kmax = 0.3 Mpc^^ is more than adequate, 
see Section 4.7. 

Since we calculate only collisional sources which are multiplied by the visibility func- 
tion, many time integrals could be restricted to 200 Mpc/c < < 500 Mpc/c. This also 
applies to the Green function integrals, where the integrand is tight-coupling suppressed 
as discussed in Section 13. However, in the numerical results presented below, we use r][^i 
as lower bound for the time integrals. The adaptive Monte Carlo algorithm (Vegas) pO|| 
automatically samples fewer points between rji^i and 200 Mpc/c than for later times. The 
upper bound of the time integrations will be denoted by rjfs. To verify that rjfg = 500 Mpc/c 
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is indeed a reasonable choice, we study the convergence of the result with increasing r/fg, 
see Section 4.7 . 

Only source terms Sx,im with small multipole moment / have to be taken into account, 
since higher first-order moments are only generated by free streaming which takes time. 
For typical values of /c, the octupole reaches its first maximum only after rj^ec^ so that 
its contribution to the integrand which is peaked around recombination is small. Higher 
moments are generated even later. If Eq. ( 3.11[ ) is used to compute the polarization modes, 
there is an additional suppression of the intensity sources Sj^im'- since the function jnm 
mixes only the two polarization modes, the intensity source contributes to polarization 



solely via the Green functions in the second line of Eq. (3.11). But as stated in Section [4.3| , 
a Green function Gmn(^)^') for multipoles coupled only by the scattering term vanishes 
if the scattering rate |k| goes to zero. Accordingly, due to the additional \k\ factor in 
front of the Green function, for small scattering rates the contribution of Sj^im sources is 
proportional to \k\^ and thus suppressed much earlier than other contributions proportional 
to \k\. Hence, sources S'/^/m with / > 2 are not only tight-coupling suppressed at early times 
but in addition do not contribute at late times. 

In the following we denote by n; the highest source-term multipole moment included. 
All sources Sx,im with / > ni are neglected. By the reasoning above we conclude that 
ni = 3 should already provide reasonable results, since Sx,3m is the highest source that 
contains the first-order quadrupole. This will be verified in detail in Section 4.7. Note that 
the source Sx,nim contains first-order solutions with I = ni + 1, i.e., for = 3 we need the 
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Figure 5: The ki integrand of {SE.22{k, ri)S'^ 22 i^' ^ v')) (cf- Eq. ( 3.16 )) integrated over k' . In the 
convolution integral over fci we have only performed the integrals over the orientation, not over the 
magnitude fci which is plotted on the abscissa, for 77 = 77' = 290Mpc/c and (from top to bottom 
in the second peak) k = 0.02Mpc"\ k = 0.04Mpc"\ k ^ 0.06 Mpc"^ and k ^ 0.08Mpc"\ The 
primordial power-spectrum is set to P'^ik) — 1/k^. The main contribution is always located around 
ki ss k. If fci is larger than fc, the integrand is suppressed by two primordial power spectra with 
large arguments. This allows us to cut the fci integration at some fci sufficiently larger than fc. 
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X 


tight-coupling solution 


numerical value 


I 


8063 ± 9 


8128 ± 10 


E 





14.6 ±2.9 


B 





0.01 ±0.02 



Table 1: Comparison of the tight-coupling sohition and the numerical result for the second-order 
power spectrum ^22 2 (^); ^ — 0.05 Mpc , -q — 200Mpc/c, primordial power-spectrum P,s>{k) = 

first order solutions up to / = 4. 

We compute all first-order quantities and Green functions in advance, before perform- 
ing the integration. However, in the integrand they are required as function of conformal 
time and wave number. Thus we tabulate the solutions, storing their values for equidis- 
tant wave-number values /cq = 0, . . . , kn,.-! = femax and Ua time values tjq = r/ini, . . . , 
rjna-i = Vfs- The time steps are chosen such that the corresponding scale factors are 
equidistant. We obtain one x Ua array for each first order quantity and one rikXnaX ria 
array for each Green function. In the further computation, the tables are interpolated 
using basis splines (B-splines). To obtain accurate results, the distance between successive 
rj and k values has to be much smaller than the typical scale of the interpolated function's 
fluctuations. 

4.6 Test in the tight-coupling regime 

In order to test our numerics we compare the numerical results for the two-point function at 
early times before recombination with the analytically known second-order tight-coupling 
solution. At early times polarization vanishes and the photon quadrupole is given by 
10Cm'i,mV^^^^_^^{ki)v^j^^^{k2)- Table ^ shows the analytical and numerical values of the 

tensor power-spectrum P22 2 (^) ^ ~ 0.05 Mpc~ at conformal time t] = 200Mpc/c. 
For this test, the primordial power-spectrum is set to P^{k) = 1/k^. As expected, the 
polarization modes are strongly suppressed and the unpolarized quadrupole has the correct 
size. The small difference between the numerical values and the tight-coupling solution 
comes from the finite scattering rate limiting the validity of the tight-coupling solution 
for which an infinite scattering rate is assumed. Evidently, in contrast to the unpolarized 
quadrupole, the E mode is strongly suppressed. Even stronger suppressed is the B mode 
since it is not coupled directly to the unsuppressed intensity quadrupole but only indirectly 
via the E mode by the free-streaming term. 

4.7 Numerical stability tests 

The numerical error of our final result is made up of a systematic error and the statistical 
error from the Monte Carlo integration. While for the latter the Vegas integration routine 
provides a reliable estimate, quantifying the systematic error requires involved analysis. 

The systematic error is controlled by five parameters: r/fg, /cmax; na, and n;. Larger 
values for any of the five parameters yield better results but require more CPU time. 
To balance both demands, we study the variation of the numerical result under variation 
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Figure 6: Dependence of the numerical result on the numerical parameters. The tensor (m = ±2) 
contribution to the value for I — 100 (top row) and for I — 1000 (bottom row) is plotted 

vs. (from left to right) the number of scale- factor steps Ua and wave-number steps nk, the number 
of multipole moments taken into account in the sources 5„, the scale factor afs from which we 
start using the free-streaming solution, and the truncation of the wave- number integrals fcmax- AH 
parameters except for the parameter used as abscissa are kept fixed at the values Ua = Uk = 100, 
= 3, rjfs = 500Mpc/c (corresponding to Ofg — 0.002064) and fcmax — 0.3Mpc~^. The central 
value of the result obtained using the default parameters is indicated by dashed lines. 



of the parameters. Values which yields a systematic error on the few percent level are 
"-a = n-k = 100, ni = 3, rjfs = 500Mpc/c and /cmax = 0.3Mpc~^. 

Figure ^ shows the stability of the result for the tensor contribution [m = ±2) to the 
angular power spectrum computed using these parameters. In the figure we identify the 
choice of r/fs to dominate the systematic error for small /, whereas for large / neglecting 
higher multipoles in the sources is the most relevant approximation. This can be easily 
understood by considering the / dependence of the spherical Bessel- function ji{k{r]Q — -q)) 
in the line-of-sight integral. For given /, the Bessel-function is small if k{riQ — rf) < I and 
oscillates for /c(r/o — rj) > /; the main contribution to the integral comes from k{rjQ — r]) ~ /. 
Since % ^ ??fs > ??, the integral is therefore dominated by contributions at k w l/rjQ. 
Accordingly, the precondition for free streaming, <C k, is fulfilled at later times for 
small / than for large /. Consequently, the computation for small / requires a large r/fg. 
With increasing I, the main contribution to the integral moves to larger k values. While 
for small k the lowest multipole moments of the first-order result are largest, with increasing 
k they are superseded by the higher moments. This can already be noticed in Fig. |6|, where 
increasing the number of multipole moments in the source terms from = 3 to = 4 
or 77-/ = 5 changes the result significantly for / = 1000, but not for / = 100. However, 
Fig. shows that for both, the tensor and the vector contribution to the C^^ B-mode 
angular power spectrum this effect does not gain importance for even larger /: the spectra 
for ni = 3 are close to the ones computed with = 4 and 5 for all I, whereas n/ = 2 is 
clearly insufficient to cover the entire range up to / ~ 2000. 
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Figure 7: Tensor (|m| = 2, upper plot) and vector contribution (|m| = 1, lower plot) to the C^^ 
spectrum, computed with different values of n;, the number of multipoles in the source terms. 



5. Result and conclusion 



We summarize our results by showing the B-mode angular power spectrum Cj^^ from the 
second-order collision term in Fig. ^ both, on a logarithmic (upper plot) and on a linear 
scale (lower plot). Similar to the weak-lensing induced contribution, the angular power 
spectrum grows rapidly with / until it reaches a maximum around / ~ 1000. At / < 200 
the tensor mode contribution (|m| = 2) is larger than the vector (|m| = 1) contribution, 
but the latter dominates in the peak region. The vector mode also shows a double peak 
structure around / 1000. 

We find that for |m| = 1 the source terms Sn are dominated by the terms containing 
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Figure 8: The C^^ angular power spectrum from second-order scattering sources (points) 
compared to the spectrum induced by primordial tensor perturbations with tensor-to-scalar ra- 
tio r = 10^^ (dashed line) and the weak-lensing signal cleaned by a factor 40 (solid line, upper 
panel only). We display separately the second-order contributions from vector (|m| = 1) and ten- 
sor (|m| = 2) perturbations, and their sum. The upper and lower plot show the same numerical 
data, first on a logarithmic, then on a linear scale. (Numerical parameters are fixed to the values 
iT-a = nk = 100, ni = 4, rjfs = 500Mpc/c, and fcmax = 0.3Mpc~^.) 



a product of the first-order electron velocity v^^^^-^ and radiation dipole A^^]^ (closely 
related to the electron velocity in tight coupling) with each other or themselves. Setting all 
other first-order quantities in Sn to zero, only a few terms remain. The source term Sj^im 
consists of only the fourth and seventh line of the right-hand side of Eq. ( 2.18 ), the source 
SE,im of the third and fourth line of the right-hand side of Eq. ( 2.19| ), and the sources 
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SB,im sjoA Sy^^im Vanish altogether as can be seen from Eq. ([23c| ) and Eq. (l2^ . With 
these approximations, we obtain the double-peaked spectrum with the first and second 
peak reaching more than 80% and 60%, respectively, of the corresponding peaks in the full 
result, Fig. |^. This indicates that the double-peak structure is related to the evolution of 
the first-order electron velocity /radiation dipole. The electron velocity oscillates over time 
until recombination and thus throughout the time interval relevant for the integration. The 
frequency depends on the wave number which in the source terms is ki and k2, respectively. 
For the vector-mode contribution the wave-vector convolution integral is dominated by 
contributions at A;i ~ A; and ^2 small or vice versa and the k integral by contributions at 
k ~ l/rjQ. Since r/o ~ 14000 Mpc/c, this relates the I values of the two peaks, I ~ 750 and 
/ ~ 1250, to the wave numbers k ~ 0.05 Mpc^^ and k ~ 0.09 Mpc~^. For both k values, 
the oscillation of the electron velocity reaches an extremum at recombination, where for 
k ~ 0.09 Mpc~^ it has undergone one full oscillation more than for k ^ 0.05 Mpc"^ The 
intermediate wave number k ~ 0.07 Mpc~^, for which the oscillation stops at recombination 
with a phase shift of only half a period, corresponds to / ~ 1000. At this / value we observe 
a local minimum of the spectrum. The tensor spectrum does not show a double-peak 
structure, since a larger range of wave numbers contributes to the convolution integral so 
that any oscillations are averaged out after integration. 

The amplitude of collision-induced B-mode polarization is several orders of magnitude 
smaller than the weak-lensing signal. At I w 200 it is comparable to the B-mode power 
generated by primordial gravitational waves with a tensor-to-scalar ratio r ~ 10~^ (shown 
as dashed line in the Figure), growing to 2 • 10"'^ at / ~ 1000. The small amplitude of 
the collisional contribution is presumably due to the fact that the sources involve only the 
cross-coupling of the first-order photon perturbations to the velocity of the baryon fiuid 
and further are mainly localized in time to the recombination era. It may therefore be 
of interest to include the reionization history into the computation to obtain a realistic 
estimate of the power spectrum at low values / ~ 10. 

The collision-induced spectrum may be compared to the spectrum induced by second- 
order vector and tensor metric perturbations ||l^. These act as sources for vector and 
tensor perturbations of the radiation intensity, which can then be converted into E-mode 
polarization by Thomson scattering through the first-order collision term, and subsequently 
to B-mode polarization by free-streaming. We find that both contributions, the present one 
and the one from metric perturbations, are similar in magnitude. If the weak-lensing signal 
could be "cleaned" completely, these second-order terms would constitute a background to 
the search for primordial gravitational waves at the level of r < 10"^ ...2- 10^^ (depending 
on which represents a challenge to CBR polarization experiments. 
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A. Summary of equations 

In this appendix we summarize further equations related to the content of the main text. 
For conventions and notation not stated exphcitly we refer to Ref. |jll|]. The equations 
are given in Fourier space conjugate to position, vectors and tensors are decomposed into 
spherical components. Phase space densities are expanded into multipoles with respect to 
the photon momentum direction. 

In the absence of first-order vector and tensor perturbations the first-order Boltzmann 
hierarchy for the energy-integrated photon phase space density multipoles (|2.3| ) is given in 
conformal Newtonian gauge by 

(A.l) 

Since B-mode polarization is not generated at first order from scalar perturbations, we may 
set A^^;^(A;) = in these equations. The Boltzmann equations for the massless neutrino 

variables A[,^^'*^(fc) are the same as the equation for A^^^'*^(fe) with the collision term on 
the right-hand side set to zero. Baryons and electrons form a tightly coupled fluid with 
first-order variables = [5 ph / Pbf^\k) , ^eVml^^-' determined by 



e,[m\ ^ 

m=— 1 



1^ + He) v^^l^ik) + zfe[-U(i)(fc) = J3. {4.«„,(fc) - A;;L(fc)} (A.2) 



with R = 3/>b/(4/)^). The cold dark matter perturbations 5^\k) = [5pc/ Pc\^^\k), v^;^^^]^{k) 
are described by the same equations with the collision term on the right-hand side set to 
zero. The system closes with the Einstein equations for the first-order metric perturbations, 
which can be put into the form 

k^D'-^\k) + 3HcD^^Hk) - 3H^A'-^\k) = Ai^Ga^ 5p^^\k), 
C;;,fn.kirm]hm2] (^^'^ + D^'^m = 87rGa^am^^^{k). (A.3) 
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The density perturbation and anisotropic stress are given by 



- ^ a(i) , a(i) 



= 15^ l^"^^A2n. + ~P^^'^:L) ■ (A.4) 

In practice, we solve the first-order equations only for mode vectors k aligned with the 
3-axis, in which case /c''^] = ik and /c'^^l = 0. In the absence of vector and tensor modes the 
m = ±l,zb2 components all vanish for this choice of k. The solution for general k needed 
in the source term of the second-order equations is obtained from the m = solution for 
k = ke^ by a rotation. In general, if Tim{k) is a spherical tensor of rank I, the relation 
between its components for k = kk and k = ke^ is 



TUk) = T,^.(te3)Z)2m(^"') = Y.^im'ike,)Y,;^"^'Ck) (A.5) 



where D^,^{R ^)m'm denotes the Wigner function for the rotation R ^ with k = Re^ 



which can be expressed in terms of the spin- weighted spherical harmonics |21|. We apply 
this to D^^\ aJ^Jq, 5\^^ etc. with / = 0, to A^„, v^^^^^ etc. with / = 1 and to aJ^,^ 

and A.^y\j^ with any /. In this case, due to the absence of first-order vector and tensor 



perturbations, only the m' = term contributes to the sum in Eq. (A.5). 

Since we do not consider polarization induced by the second-order vector and tensor 
metric perturbations, themselves induced by the first-order scalar modes, we do not need 
to solve the second-order Einstein equations. The only second-order quantity other than 
the photon perturbations that we must solve for is the electron velocity for which we use 
the truncated equation ( ^.211) as explained in the main text. For completeness we provide 
here the complete equations for the second-order variables of the baryon-electron fluid. We 
define the fluid variables through the energy-momentum tensor in the local inertial frame 
at rest and aligned with the general coordinate system 

Tab = [eAT[(iBY T^y = {p + p)uAUB -PVAS + ^ab 

■ g{r],x,p) (A.6) 



(2^)3^'" E 

Here p"^ denotes the physical (not co-moving) particle momentum in the local frame and 
g{ri,x,p) the phase-space distribution, for which we take a local Maxwell-Boltzmann dis- 
tribution. The four- velocity is parameterized in the form = (1/^1 — i?^, v). We then 
calculate the conformal-time derivative of the expression in the second line using the results 
of Ref. Q- Together with 

Too = P+-p[v^^^? + ..., To, = -pv' + ..., (A.7) 
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valid to second order in perturbations, we obtain the desired equations of the density and 
velocity perturbation of the baryon-electron fluid to second order: 



m=— 1 



■n=-l ^ ^ 



+ (-l)™^c<[L](fci)<[lH(^2)-6Z)«(fci)i)«(fc2) + 3I)«(fci)5W(fc2 



m=~l 



m=— 1 

(1) 



m'=-l 



.(1) 



(1) 



,(1) 



,["12] 



(A.8) 



(fci)(4-SL](fc2)-A(;L(fc2)) 



(A.9) 



The left-hand sides of these equations are obtained from the I = and I = 1 moments of the 
Boltzmann equation for massive particles given in |11]. The right-hand sides follow from 
the collision term for photons using the tight coupling of baryons and electrons through 
Coulomb scattering and energy-momentum conservation in photon-electron Compton scat- 
tering. Eqs. ( |A.8| ) and (A^) agree with the corresponding results in Ref. [^]. 



Finally, the coupling coefficients Cmim2 and Dmim2 that appear in the Boltzmann 
equations read 



a 



m±l,m 



a 



m±l,m 



c- 



^J{l + l±m){l + 2±m) 
V2{21 + 3) 



21 + 3 
V2{21 - 1) 



21 - 1 
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0,/ _ J2{l + l±m){lTm) 

^m±l,m r 



^(^ + 1) 

Q I 2m 
B. Derivation of the line-of-sight solution 

The line-of-sight sohition of the first-order Boltzmann equations have been derived in 
Ref. [^4| (see also Ref. fl^). We provide here a generalization which allows for more 
general source terms including higher multipole moments and polarization modes. 
Our aim is to solve the equation 

+ kCnm^m — ~\f^\^n + Pn 

for k = ke^. First we solve the homogeneous differential equation without the source 
term p„. This can be done easily after reversing the multipole decomposition. Then the 
homogeneous equation takes the simple form: 

Aab + in ■ k Aab = -\k\Aab, (B.2) 

where ab = ++, H — , — h, are the polarization indices. This equation is solved by 

Aabiv) = e--'=(^-^')-''('''^') AabW), (B.3) 

where K{r],r]') is the integral over \k\ from r]' to r]. If the first argument of k is the 
conformal time today ijo, we omit this argument in the following and in the main text, and 
write K{r]) = K{r]Q,r]). Next, we decompose this solution into multipoles. For this purpose, 
we expand the exponential in spherical harmonics and Aq{,(V) ™ spin-weighted spherical 
harmonics, 

h 



h,rn 




and apply on both sides i^\J'^^^ j d^^imin) to obtain 

^ Zi2 + i 

The spin s is zero for ab = ++, and ±2 for ab = ib=p. Finally, instead of the photon 

helicity state densities, we use the Stokes parameter distributions 

^I,lm = ^(A++,;m + A ^im), (B.6) 

AE^im = 2 (A+-,«m + A_+_;m), (B.7) 

Ab,Zto = -(A_|__^;m — A_+_;m)- (B.8) 
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On the right-hand side of Eq. (|B.7] ) and Eq. ( [B.8D there is a sum and a difference, respec- 
tively, of Clebsch-Gordan coefficients with opposite spin. They are related by 

(-Vo')-<-')'-"-'^(,:'o':)- («-^' 

Depending on the parity oi I — h — I2 we obtain either a sum or a difference of the H — 
and — + components. For odd parity, this leads to a mixing of E and B modes. The 
resulting mode coupling can be written in a compact form as matrix Hxx' defined by 
Hxx'ieven) = 5xx' and Hxx'iodd) = 5xiSx'i + i^XE^x'B - iSxB^x'E- Then 



h,h X'=I,E,B 

m m j \Fx F: 



X' 



H*xx'{l-h-h)^x'MM), (B.IO) 



where Fj = 0, Fe = Fb = —2. A more detailed explanation of the effect which introduces 
the matrix H and the properties of spin-weighted spherical harmonics can be found in 
Ref. |11|. With the solution of the homogeneous equation at hand, we can immediately 
write down the solution of the inhomogeneous equation: 



ii,i2 X'=I,E,B 



H*xx'{l-h-l2)px',i,M)- (B.ll) 



\ m m I \ Fx Fx' 
Using the multi-index notation from the main text, we can define the function 

^(,) i^ln + I)i2h + ^) ^^^^^^H*x„xJln - h - Im) 



Jnm\ 

' h 



mn m„J \Fx„ Fx„ 

and P{vc,im) = to obtain the compact equation 

•no 



(B.12) 



An(r?o)= / dr^e-''^'^^jnm{Km-r]))pm{r]), (B.13) 



which holds for X„ = I, E, B. 
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